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, Abstract 

Several concepts on the measure of observability, reachability, and 
robustness are defined and illustrated for both linear and nonlinear 
I control systems. Defined by using computational dynamic optimiza- 

1/-^ ■ tion, these concepts are applicable to a wide spectrum of problems. 

Some questions addressed include the observability based on user- 
I information, the determination of strong observability vs. weak ob- 

' servability, partial observability of complex systems, the computation 

■ of L^-gain for nonlinear control systems, and the measure of reachabil- 

^\ , ity in the presence of state constraints. Examples on dynamic systems 

defined by both ordinary and partial differential equations are shown. 



1 Introduction 



Control systems are analyzed and characterized by using fundamental con- 
cepts such as observability, reachability, and input-to-output gain [H [21 13] . 
These concepts have a vast volume of literature. For nonlinear systems, the 
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challenge is to define the concepts so that they are characteristic and fun- 
damental to control systems and, meanwhile, they are practically verifiable. 
In this paper, the goal is to use dynamic optimization to define quantitative 
measures of control system properties. Moreover, computational methods of 
dynamic optimization provide practical tools to numerically implement these 
concepts in applications. 

In Section [21 the ambiguity in estimation is defined as a measure of ob- 
servability. This quantity can be numerically computed by solving a dynamic 
optimization. An example is shown in which some systems observable in tra- 
ditional sense are not practically observable because of their poor value of 
ambiguity in estimation. In other words, we can quantitatively tell strongly 
observable from weakly observable. Another feature of this concept is the 
capability of taking into account non-sensor information or user knowledge 
of systems, in addition to the output. For instance, an example is shown in 
which the system is unobservable under a traditional definition. It turns out 
that the system is strongly observable if we know the control input has a 
bounded variation, but without using its accurate upper bound. Moreover, 
this concept can be used to measure partial observability of complex systems, 
including the observability of a function of the states and the observability 
of unknown parameters in a model. 

In Section [3l computational methods for the L^-gain of control systems are 
introduced. The assumption is that the space of input has finite dimension. 
Then, the L^-gain can be computed using dynamic optimization. In addition, 
a method of approximating L^-gain is also introduced, which is based on 
the eigenvalues of covariance matrices. The methods are exemplified by a 
nonlinear model of atomic force microscope. 

In Section HI we define the concepts of ambiguity in control and control 
cost. These definitions take into account the control input as well as systems' 
constraints. For instance, the concept can be used to quantitatively measure 
the reachability of nonlinear systems under the constraint that the states 
must stay in a given region of safety. As an example, heat equation with 
boundary control is studied. 
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2 Observability 



Consider a general control system 

G M G 3?"" /i G 3?"^ 

(1) 



X = f{t, x,u,iJ,), X G 3fJ"", M G 3?"" /i G JJ"^ 
?/ = x,u,fi), ye S?"'^ 
2; = e(t, X, M, /i), 2; G 3?"^ 



in which |/ is the output, z is the variable to be estimated, which is either the 
state a; or a function of x in the case of partial observability for large scale 
systems. The system state is x, u is the control input, /x is the parameter or 
model uncertainty. In ([T]), C is a general formulation of constraints. Some 
examples of constraints include, but not limited to, 

E{x{tQ),x{tf)) < 0, end point condition 
s(x,m) < 0, state-control constraints 

s{x{ti)) = 0, known event at time ti 

l^min < Ai < fJ'max modcl Uncertainties 

s{x,fi) = 0, DAE (differential-algebraic equations, // is a variable) 

control input with bounded variation 
(non-sensor information) 



Variationfw) < 



These constraints represent known information about the system in addition 
to the measured output y. This general form of constraints makes it pos- 
sible to take into account non-sensor information, or user knowledge about 
the system, in the estimation process. For instance, some state variables are 
known to be nonnegative; or a control input has bounded variation; or an 
event is known to happen at certain moment. All these are valuable infor- 
mation that can be used for the estimation of z. The goal of this section is 
to define a measure for the observability of z using the observation data of y 
as well as the constraints and the control system model. 



2.1 Definition 

We assume that variables along trajectories are associated with metrics. For 
instance, y = h{t,x{t),u{t), fi), as a function of t, has or norm; z = 
e{t,x(t),u{t), fj,) can be measured by its function norm, or by the norm of 
its initial value e(to, ^(to) , u(to) , fi) . A metric used for z is denoted by || • ||z; 
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and 1 1 ■ I |y represents the metric for y — h{t, x, u, fi). The following definition 
is applicable to systems with general metrics, including U' and L°°. Unless 
otherwise specified, a norm | |a|| for a e 3?*^ is defined by 




which equals its essential supremum value. In this paper, a triple {x(t),u(t), fi) 
represents a trajectory of (1) satisfying the difi'crcntial equations as well as 
the constraints. Given a positive number e > and a nominal, or true, 
trajectory {x(t),u(t), fi). Define 



S = {{x{t),u{t),ft)\ \\h{t,x{t),u{t),ft) - h{t,x{t),u{t),fj,)\\Y < e} (2) 



The number e is used as an output error bound. If h{x{t),u{t), p,) stays in 
the e neighborhood of the nominal output h{t, x{t),u{t), n), then we consider 
the trajectory {x{t),u(t), fi) not distinguishable from the nominal one using 
output measurement. In this case, any trajectory in £ can be picked by an es- 
timation algorithm as an approximation of the true trajectory {x{t),u{t), fi). 
For this reason, a trajectory in £ is called an estimation of {x{t),u{t), n). 
Similarly, z — e{t, x, u, fx) is an estimation of 2; = e(i, x, /i). 

Definition 1 Given a trajectory {x{t),u{t), fi), t e [to,^i]- Let e > be the 
output error bound. Then the number Poi^) is defined as follows 




ial + ... + aiy/' 



For any function h{t), t G [to,tf], its L^-norm is defined by 




The infinity norm is defined by 




Po{e) 
subject to 



max 

(i(t),u(t),/i) 



e{t,x{t),u{t),ii) - e{t,x{t),u{t),ii)\\z 



\\h{x{t),u{t),fi) - h{t,x{t),u{t),ii)\\Y < e 

X = f{t,X,U,fl), 

(£(•),«(•), A) eC 



(3) 
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The number Po{^) is called the ambiguity in the estimation of z along the 
trajectory {x(t),u{t), ft) . 

Let U be an open set in {x,u, fi)- space and [to,ti] be a time interval. 
Then the largest value of ambiguity along all trajectories in U is called the 
ambiguity in the estimation of z in the region U . 

Remarks 

1. The ratio po{t)/e measures the sensitivity of estimation to the noise in 
y. A small sensitivity value implies strong observability of z in the presence 
of sensor noise. 

2. The ratio po{t)/e is closely related to the observability gramian. Con- 
sider a linear system 

X = Ax, y = Cx 

Suppose z = X and suppose 1 1 ■ | |y is the L^-norm. Let P be the observability 
gramian [1], [1], then 

\\y\\l = xlPxo (4) 

Given ||?/||y = e, po equals the maximum value of ||xo|| satisfying (jl]). In this 
case, Xq is an eigenvector of P associated to the smallest eigenvalue Amm- We 
have 

= ^minpl (5) 

Therefore, the ratio po{e)'^/e^ equals the reciprocal of the smallest eigenvalue 
of the observability gramian. For nonlinear systems, one can use observability 
gramian to approximate po(e)/e. An advantage of this approach is that 
the gramian can be computed empirically without solving the optimization 
problem ([3]). Details on empirical computational algorithms for the gramian 
of nonlinear systems can be found in [1] and [5]. 

3. Given a fixed number po > 0, from (jlj) the least sensitive initial state 
defined in ([5]), i.e. the eigenvector of length po associated to Xmin, can be 
found by the following optimization 

arg min 1 |y = Xq Pxq 

\\xo\\x=Po 

Extending this idea to nonlinear systems, the least observable direction in 
initial states can be defined as follows: given po > 0, let e be the minimum 
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value from the following problem of minimization 

e = ,,,nim \\h{t,x{t),u{t),fi) - h{t,x{t),u{t),fi)\\z 

subject to 

ll^^o - a^olU = Po 
X = f{t,x,u, p,), 
{x{-),u{-),fi) G C 

The resulting initial xq represents the least observable state on the sphere 
lla^o — a^olU = Po and the ratio po/e measures the unobservability of initial 
states. This definition is a reverse process of Definition [H However, in some 
cases it is easier to handle the constraint | |a;o — Xo| |x = Po than the inequality 
of h(t, X, u, /i) in ([2]). 

4. The metric for output in Definition [1] can be a vector valued function 
which is bounded by a vector e. This flexibility is useful for systems using 
different types of sensors with different accuracy. 

5. Definition [T] is independent of estimation methods. It characterizes 
a fundamental attribute of the system itself, not the accuracy of a specific 
estimation method. In the following, we compare Definition [1] to traditional 
definitions of observability. It is shown that simple linear systems observable 
in the traditional sense might be weakly observable or practically unobserv- 
able under Definition [H and, on the other hand, some systems not observable 
under traditional definitions are practically observable with a small ambigu- 
ity in estimation, o 



2.2 Computational dynamic optimization 

The problem defined by ([3]) is a dynamic optimization. To apply Definition 
[H this problem must be solved. Obviously, an analytic solution to ([3]) is 
very difficult to derive, if not impossible, especially in the case of nonlinear 
systems. However, there exist numerical approaches that can be used to 
find its approximate solution. For instance, various numerical methods are 
discussed in detail in [6J, [7J, and [8j. Surveys on numerical methods for 
solving nonlinear optimal control problems can be found in [9l [10]. The 
computational algorithm used in this paper is from a family of approaches 
called direct method |TT[ [121 HSl [HI [IS] . The essential idea of this method is 
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to discretize the optimal control problem and then solve the resulting finite- 
dimensional optimization problem. The simplicity of direct methods makes 
it an ideal tool for a wide variety of applications of dynamic optimization 
with constraints, including ([3]) in Definition [H 

More specifically, all simulations in this paper use a pseudospectral opti- 
mal control method. In this approach, a set of nodes is selected using either 
the zeros or the critical points of orthogonal polynomials, in our case the 
Legendre-Guass-Lobatto nodes. Then, the problem of dynamic optimiza- 
tion is discretized at the nodes to result in a nonlinear programming, which 
is solved using sequential quadratic programming. Details are referred to 
fr2i [TU [12] . In some of the following examples, dynamic optimizations are 
solved using the software package DIDO [TB]. 

A frustration in nonlinear programming is the difficulty of finding global 
optimal solutions within a given domain. This is no exception in this paper. 
In all examples, a variety of initial guesses are used to gain a comfortable 
level of confidence that the result is not a local optimal solution. However, 
for all examples of nonlinear systems in this paper, the computation cannot 
guarantee global maximum value for Nevertheless, in the case that a 
result is not the global maximum value, it still provides a lower bound of the 
ambiguity value po- 

2.3 Examples 

For the rest of this section, we illustrate Definition [1] using several examples. 
In the first example, it is shown that the traditional concept of observability 
is ineffective for systems with large dimensions. It justifies the necessity of a 
quantitative definition of observability, such as Definition [TJ 

Example. Consider the following linear system 

Xi = X2 
X2 = Xs 



y = xi 

Under a traditional definition of observability, this system is perfectly observ- 
able for any choice of n, i.e. given an output history y = Xi{t), it corresponds 
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to a unique initial state Xq. However, if Definition [T] is applied to measure the 
observability, it is a completely different story when the dimension is high. 

Suppose the goal is to estimate xq. We can define z = x{t). Definition 
[1] is applicable with arbitrary metrics. To measure the observability of the 
initial state, we can use the norm of x(0) as the metric for z, i.e. 

lkWllz=||a:(0)|| 

For this example, the output accuracy is measured by L°°-norm, 

umy 

Let us assume that the true initial state is 



max \v(t)\ 



Xo 







Let the output error bound be small, e = IQ-^. So, we assume very accurate 
observation data. The time interval is [0, 15]. Problem ([3]) has the following 
form 

Po(e) = max ||x(0) — a;(0)|| 

X 

subject to 

||Xi(t) -Xi(t)||y < e 

X = f{x) 

It is solved to compute Po(e)- Table [1] lists the result for n = 2, 3, ■ ■ ■ , 9. 



(7) 



n 


2 


3 


4 


5 


Po(e) 


4.70 X 10-6 


2.67 X 10-5 


1.53 X 10-^ 


8.89 X 10-^ 


e 


10-6 


10-6 


10-6 


10-6 


n 


6 


7 


8 


9 


Po(e) 


5.20 X 10-3 


3.01 X 10-2 


1.75 X 10-1 


1.02 


e 


10-6 


10-6 


10-6 


10-6 



Table 1: Observability 



From the table, when n = 2 the ambiguity in the estimation of xq is as 
small as 4.70 x 10-6. So the system is strongly observable. Equivalently, if 
the observation data has absolute error less than e, then the worst possible 
estimation of xq has an error at the scale of 10-6. This conclusion agrees 
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Figure 1: Estimation error (n = 9) 



with the traditional theory of observabihty. However, when the dimension is 
increased, the observabihty ambiguity increases too; thus the system becomes 
less observable. At n = 8, the observability ambiguity is as big as 0.175, or 
the worst error of estimation is 17.5% relative to the true xq- When n = 9, 
the observability ambiguity is 1.02. In this case, the worse relative error in 
estimation is more than 100%! Thus, the system is practically unobservable, 
although it is perfectly observable under a traditional definition. Figure [1] 
shows why this system is practically unobservable. The continuous curves 
represent the true trajectory and its output for n = 9; the dotted curves are 
the estimation. The outputs of both trajectories agree to each other very well 
(Figure on top), but the initial states (only xg is plotted) are significantly 
different. 

As shown in Figure [T|, while the estimation of the initial state is inaccu- 
rate, the estimation is very close to the true value at the final time tf. To 
see the observability of the final state, let us use a different metric for z, 

Mmz = Mtf)\\ (8) 

If we consider t = tf a.s the current time moment, then this metric is used 
to measure the detectability of the current system state, rather than the 
observability of the initial value xq- To compute the ambiguity under the 
new metric, we solve the problem defined in ([7]) except that the cost function 
is replaced by the metric ([HD. For the case of tf = 10 and e = 10~^, the 




2 4 6 8 10 12 14 16 

time 
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ambiguity in the estimation of x{tf) equals 2.7328 x 10~^. Therefore, the 
system is accurately detectable. 

To summarize, this example shows a set of linear systems that are observ- 
able under conventional definition. However, as the dimension is increased, 
the systems become practically unobservable in the sense that an output 
trajectory cannot accurately determine the state trajectory. Meanwhile, the 
detectability of the system is not changed with the dimension. Definition [T] 
is used here to treat both observability and detectability in the same frame- 
work, quantitatively. <) 

A concept is useful only if it is verifiable for a wide spectrum of systems 
and applications. An advantage of Definition [T] is that the dynamic optimiza- 
tion ([3]) can be numerically solved for various types of applications. In the 
following, we illustrate the usefulness of the ambiguity in estimation using 
two examples. One is a networked cooperative control system; and the other 
one is parameter identification for nonlinear systems. 

Example (Partial observability of cooperative and networked systems). In 
this example, it is shown that an unobservable system under traditional con- 
trol theory can be practically observable by employing user knowledge about 
the system, such as an approximate upper bound of control input. Con- 
sider a networked control system showing in Figure [2l Suppose it consists of 
an unknown number of vehicles. Due to the large number of subsystems, it 
could be either impossible or unnecessary to process or collect all information 
about the entire system. A practical approach is to find partial observabil- 
ity with local sensor information only. In this example, we assume that the 
cooperative relationships in the system is unknown except that we know Ve- 
hicle 2 follows Vehicle 1; Vehicle 3 follows both Vehicle 1 and 2 as shown by 
the arrows in Figure [21 The dashed lines in the figure represent unknown 
cooperative relationships. The question to be answered is the observability 
of Vehicle 1 if the locations of Vehicle 2 and 3 can be measured. 

Suppose each vehicle can be treated as a point mass with a linear dynam- 
ics 



Assume that the control input of vehicles 2 and 3 are defined as follows 



Xil = Xi2 Vil = yi2 
Xi2 = Ui yi2 = Vi 



U2 = ai(x2i — xu — di) + a2(x22 — 2:12) 
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Figure 2: Cooperative networked system 



where di is the distance of separation. The control in the ?/-direction is the 
same. So, Vehicle 2 follows Vehicle 1, Vehicle 3 follows the average position 
of Vehicle 1 and 2. Suppose we can measure the positions of Vehicles 2 and 
3. 

output = [ X21 1/21 xsi ysi J (9) 

The question to be answered is the observability of the location and velocity 
of Vehicle 1, i.e. Xn, yu, Xi2, yi2- We would like to emphasize that the 
control input of Vehicle 1 is unknown because its input is determined by its 
cooperative relationships with other vehicles or agents, which is not given. 
Therefore, in traditional control theory. Vehicle 1 is unobservable. 

To make Vehicle 1 practically observable with limited local measurement, 
we assume that the input of Vehicle 1 has bounded variation with an upper 
bound Vmax- This is to say that the vehicles are not supposed to make high 
frequency zigzag movement, or the control does not have chattering phe- 
nomenon. However, discontinuity in control, such as bang-bang, is allowed. 
In the following, we measure the observability of Vehicle 1 along a trajec- 
tory defined by the parameters in Table [21 The control input of the nominal 
trajectory is 

. {tf-to)t 
Ui = sm — 

TT 

which is unknown to the observer. To measure the ambiguity in the esti- 
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di 


d2 




02 


&1 


^2 




1-^21' -^22/ 


1-^31' •^32/' 





20 


-2 


-2 


-1 


-2 


-3 


-7 


(0,4) 


(rfi,4) 


(^2,4) 



Table 2: Parameters of Nominal Trajectory 



mation of vehicles, we assume the output error bound of ([9]) is e = 10~^. 
For the unknown ui, we assume a bounded variation of less than or equal to 
Vmax = 3.0, which is 50% higher than the true variation. The metric for each 
output variable is the L°°-norm. The metric for the location and velocity 
of Vehicle 1 is the L^-norm. The ambiguity in the estimation of each state 
variable is computed by solving a problem of dynamic optimization. Using 
the estimation of Xn as an example, we have 

Po(e) = max pii(t) - Xii(t)||i2 

{x,ui) 

subject to 

\\x2i{t) - a;2i(t)||Loo < ei 

\\X3iit) - X3i{t)\\Lo^ < 62 

X = fix) 

where V{u]) is the total variation. In computation, this constraint is dis- 
cretized at a set of node points to < ti < ■ ■ ■ < = ff so that 

N 

J2 \Ul{ti) - Ul{ti-.l)\ < Vmax 
k=l 

An interesting point in the formulation ffTOj) is that the outputs for the two 
vehicles, i.e. £21 and X31, have different error bounds, ei and €2- The metric 
for the outputs is a vector valued function. This flexibility of using different e 
value for multiple outputs is advantageous for systems with multiple sensors 
of different qualities that measure various states with different accuracy. 

The computed result is shown in Table [31 The small relative ambiguity 
value shows that the location and velocity of Vehicle 1 are practically ob- 
servable given the measurement of the positions of Vehicle 2 and 3, without 
using any information about the rest of the networked system and without 
knowing the input of Vehicle 1. The worst estimation of xu and xu is shown 
in Figure [3|, which has good accuracy. <) 

In the following example, the concept of ambiguity in estimation is ap- 
plied to the Laub-Loomis model [T7] with unknown parameters, a nonlinear 
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^max 


e 




Pxii/\\Xu\\l2 






3 


10-^ 


1.2257 


2.8 X 10^^ 


0.5901 


1.16 X 10'2 



Table 3: Observability of Vehicle 1 
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Figure 3: The worst estimation of the position and velocity of Vehicle 1 
system of oscillating biochemical network. 

Example (Parameter identification) In the study of biochemical networks, 
it was proposed that interacting proteins could account for the spontaneous 
oscillations in adenylyl cyclase activity that was observed in homogeneous 
populations of dictyostelium cells. While a set of terminologies such as 3', 5'- 
cycle monophosphate (cAMP) and adenylate cyclase (ACA) are involved in 
the problem, we focus on the state space in which a set of seven nonlinear 
differential equations are used as the model |T7j. 

Xi = kiXj — k2XiX2 
X3 = k^XY - kGX2X3 

Xi = kj — kgx^x^ (11) 

3^5 = kgXi — ki()X4.Xc^ 

xe = kiixi - kuXQ 
^7 = kisXe — kuXj 

In a robustness study JTE\, it was shown that a small variation in the model 
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parameters can effectively destroy the required oscillatory dynamics. As 
a related question, it becomes interesting to investigate the possibility of 
estimating the parameters in the system, ki, k2, ■ ■ ■, kr. To exemplify the 
idea, we assume that xi, the value of CAC, is measurable, i.e. 

y = xi 

We also assume that the initial states in experimentation is known. Suppose 
the unknown parameters are ki, k^, and kio; and suppose the other param- 
eters are known. The goal is to use the measured data of y to estimate the 
unknown parameters , i.e. 

z = ki ka kio 

Using Definition [1] we can quantitatively determine the observability of the 
unknown parameters. Along a nominal trajectory {x*{t),k*), the ambiguity 
can be computed by solving the following special form of 

pI = max (k, - klf + {k, - klf + {k,o - kl,f 

subject to 

\\x,{t) - ximh < 

X = fit, X, ki, kQ, kio), other parameters equal nominal value 
x(to) = x*{to) 

In the simulation, a nominal trajectory is generated using the following pa- 
rameter value and initial condition 

k^ = 2.0, k2 = 0.9, ks = 2.5, k^ = 1.5, k^ = 0.6, /cg = 0.8, k-r = 1.0, 

fcg = 1.3, kg = 0.3, kio = 0.8, kn = 0.7, k^ = 4.9, kis = 23.0, k^ = 4.5, 



x(0) 



1.9675 1.2822 0.6594 1.1967 0.6712 0.2711 1.3428 



The output error bound is being set at e = 10 ^. The computation reveals 
that the ambiguity in the estimation of z = ki kg kiQ is 

Po = 2.38 X 10-2 

Given the nominal value of the parameters, the relative ambiguity in estima- 
tion is about 1%. So, the parameters are strongly observable. In fact, the 
worst estimation of the parameters is 

ki = 2.0150, ke = 0.8082, kw = 0.7836. 

The trajectory generated by the worst parameter estimation is shown in 
Figure m o 
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Figure 4: The trajectory of worst estimation: curve - true trajectory; star - 
estimation) 

3 Input-to-output gain 

L^-gain is a tool of analysis widely used by control engineers to quantitatively 
measure the sensitivity and robustness of systems. Consider 

i = /(t,a;,w,/i) .^2) 

where x G 3?""' is the state variable, w G 3?"™ is the input that represents the 
disturbance, fi G 3?"'' is the system uncertainty or a parameter, z G 3fJ"^ is 
the performance. In the following, the L^-norm of a vector valued function 
is denoted by || ■ ||lp, for instance 




Given a fixed time interval [to;^i] and a > 0. Suppose the input w is a. 
function in space for some 1 < p < oo such that w{t) is bounded by a, 
i.e. 

||w(t)||LP < a 

Let x*{t) be a nominal trajectory with x*(to) = Xq, w*{t) = 0. Fix the initial 
value x(to) = Xq. Suppose the system uncertainty is bounded, /imm < < 
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f^max- Then the L^-gam from w to z along x*(t) is defined as follows 



7(cr) 



e(t, X, w, jj) - e(t, X*, 0, /i) | |lp 




max 

I < < fJ-n 



max 



a 



(13) 



Remark 6. Without parameter, the maximum value of | |e(t, x, w)—e{t, x*,0)\\lp 
is the ambiguity in the estimation of 2; = e{t,x,w). More specifically, con- 
sider the ambiguity in the estimation of z under the observation of w with an 
error bound a. Then the L^-gain 7(0") gain equals the ratio of the ambiguity 
and a. 

3.1 Computation and example 

The input-to-output gain can be computed by solving the problem ( IT3ll . In 
Section [2], the output function y is smooth. It can be numerically approx- 
imated in a finite dimensional space, such as interpolation using a finite 
number of nodes. Similarly, one has to work on a finite dimensional space 
of w to carry out the computation. So, for the purpose of computation, we 
only discuss the L^-gain in a finite dimensional space of w, denoted by U, 
rather than the infinite dimensional space of arbitrary integrable functions. 
The space U can be defined by the frequency bandwidth, or the order of 
polynomials, or some other spaces used for the approximation of the input. 
Then, ( ITSj) is reformulated as follows 



More specifically, given a positive number a > define 

J(a;(-),w(-),/i) = \\e{t,x{-),w - e(t, 0, |lp 

Then the following dynamic optimization determines the L^-gain over the 
space U. 




max 

W EU,\\w\\lp < (7 



e{t,x,w,^) - e(t,x*,0,fi)\\LP 



a 



(14) 
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Figure 5: Atomic force microscope 



Dynamic optimization for L^-gain 



p = max J 

{x,w,fi) 

subject to 

w{t) eU, \\w{t)\\LP < a (^i^^ 

X = f{t,X,W,IJ,), 
x{to) = Xo 

Pmin ^ p ^ pvaax 



The L^-gain from w ElA to z is 7w(o") = — • o 

o 

Example (L^-gain in the presence of system uncertainty). Atomic force 
microscope (AFM) invented two decades ago is used to probe surfaces at the 
atomic level with good accuracy. This type of equipment is also used as nano- 
manipulation tools to handle particles at nano-scale [IHl EOl El] • Illustrated 
in Figure [5l the system consists of a microcantilever with a sharp tip at one 
end. The vibration of the cantilever is measured by an optical sensor. The 
topographic images of surfaces can be taken by measuring the cantilever's 
dynamic behavior which is determined by the interacting force of the tip with 
the sample. 

The dynamics of the vibrating tip can be modeled as a second order 
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system [2T] 



Xi = X2 

X2 = —oj^x\ — 2^ijjX2 + h(x\, S) + + (16) 

Z = X\ 

where x\ is the position of the cantilever tip at the scale of nanometers, 
xi is its velocity, uj is the natural frequency of the cantilever, and ^ is the 
damping coefficient. In this system, u{i) is the control input, w{f) is the 
actuator disturbance which is unknown. The function is the tip- 

sample interaction force in which 8 is the separation between the equilibrium 
of x\ and the sample surface. It is a system uncertainty. We adopt the 
following model for h{x\, 5) [20] . 

h{xi,d) = -— — + 



(S + x^y {5 + x,f 

At nano-scale, system uncertainly and performance robustness are critical 
issues in control design because a seemingly small noise or uncertainty may 
have significant impact on the performance. In the following, we assume that 
the value of 6 and w{t) are unknown. The goal is to compute the L^-gain 
from the actuator disturbance w to the performance z = Xi for 6 in the entire 
interval [6min, Smax]- 

The simulations are based on the following set of parameter value 

uj = 1.0, ^ = 0.02, ai = 0.1481, 03 = 3.6 x 10"^ 

The nominal control input is u{t) = 1 and the time interval is [0, 7], which is 
long enough to cover one period of oscillation. The bounds of S are 

^min 0.8, SfYiax 1-2 

We use L^-norm as the metric for the actuator disturbance force w. Let w 
be an arbitrary function in a two- frequency space yVki,k2 defined as follows 

( , , 27!'ki , , 2Trkj A 

w = Y.[Ai cos(- '-t) + sin(- '-t) 

i=i \ tf-to - to J 

Let cr = 0.03. Then the L^-gain is computed by solving the dynamic opti- 
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Figure 6: L^-gain 



mization defined in (|T5l) . More specifically, 

p = max ||a;i(t) — x*(t)||i2 

{x,u,S) 

subject to 

\\w{t)\\L2 < a, w{t) e Wk,M 
X = f{t,x,u,6), 
x{to) = Xq 

The L^-gain equals — . It is computed for spaces with various frequencies, 
cr 

Wo,i, W2,3, W4,5, Wqj, and Ws,9- The result is shown in Figure El The L^- 
gain for frequencies and 1 is 2.5707. When the frequencies are increased, 
the gain decreases, o 



3.2 An alternative algorithm for L^-gain 

Solving ( JT5l) becomes increasingly difficult for high frequencies. The reason 
is that, for the computational purpose, the problem of dynamic optimization 
is always discretized at finite number of nodes in time. For higher frequen- 
cies, the number of nodes must be increased. As a result, the dimension of 
optimization variables is increased as well. Developing efficient methods of 
computation for inputs of high frequencies requires further research. How- 
ever, in the case of L^-gain, there exists an alternate approach without solving 



19 



dynamic optimization. Inspired by Remark 6 and the observability gramian 
in [1], the L^-gain can be approximated using the following matrix. 

Suppose the space of input, W, is finite dimensional with a basis, wi, W2, 
■ ■ ■, Wm- Let 0" > be a constant number. For the input ±awi, the trajectory 
of 

X = f{t,x, ±awi) 
x{to) = Xq 

is denoted by Define 



2Az' = e{t, awi{t)) - e(t, x'-{t), -awi{t)) 



Now, define 



GT, =< Wi,Wj > 



tf-t 



f —iQ -'to 



Gl^=<I\z\Az^ >=- 

If - to J to 



Wi{t)'^Wj{t)dt 
' Az\t)^Az\t)dt 



Denote the matrices = (G'^)"j=i and = (G^^.)I^i=i- Given any 



w 



i=l 



satisfying 

< W,W >= cr^ 

Then Az generated by ±w is approximately 



(17) 



(18) 



2Az = e{t, x'^{t),w{t)) - e{t, x'{t), -w{t)) ^ 2 ^ 



i=l 



Therefore, 



1 1""^! 



1 12 


'^1 ' ' ' '^m 


G"" 






1 12 ~ 

IIl2 ~ 


0-1- ■ ■ 0"m 


G' 


ai ■ ■ 





T 
T 



Therefore, the L^-gain square is approximately the solution of the following 
optimization 

— max a G a 

a 

subject to 

a^G^^a = 
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CLl 0,2 ■ ■ ■ CLy, 



where a 
tion is 

for some scalar A. At this point, 

= Xa^ 



Therefore, 

lu{(^? = „ max 



T 

. A necessary condition for the optimal solu- 



\\Az 


|2 

Il2 


\w\ 


2 

L2 



a2 



= A 

On the other hand, A is an eigenvalue of {G^)~^G^. So, the L^-gain is 
approximately the square root of the largest eigenvalue. 
To summarize, given a system 

X = f{t,x,w) 
x{to) = Xo 

and a space of input functions U with basis wi, W2, ■ ■ ■, Wm- Given a > 0, 
compute Az^ in f|T7|) . Compute the matrices G^ and G^ in f|T8|) . Then 
the L^-gain is approximately y/Xmax, where Xmax is the largest eigenvalue of 



Example. Consider the model of AFM defined in flTBl) . Assume that the value 
of 5 = 1.0 is known. The other parameters are the same as in the previous 
example. The approximate L^-gain is computed using the matrix approach. 
The result is shown in Figure [71 This method is straightforward in computa- 
tion because no optimization is required. However, the approximation does 
not take into full account the nonlinear dynamics. Comparing to the gain 
using f|T5l) , the errors of the gain computed using covariance matrix is around 
18 ~ 20% in Wo,i, Wa.s, and ^4,5- 
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Figure 7: Approximate L^-gain 



4 Reachability 

Dynamic optimization can be applied to quantitatively measure reachability. 
Consider a control system 

x = fit,x,u) (19) 

where x G 3?""" and u G 3?"". Suppose is a norm in 3?"^. We suppose 

that the state and control are subject to constraint 

Definition 2 Given xq and Xi in 3?""= . Define 

pc{xo,Xif = mill ||a;(ti) - Xi\\x 

(x,u) 

subject to 

X = f{x, u) 

(20) 

x(to) = Xq 

ix{-),u{-))EC 

The number Pc{xo,xi) is called the ambiguity in control. 

Let Dq,Di C ^^"^ be subsets in state space. The ambiguity in control over 
the region Dq x Di is defined by the following max-min problem. 

Pc = max _ pc{xo,Xi) 

(xo,xi)eDoxDi 
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In this definition, ti is either fixed or free in a time interval [to, T]. In the 
following discussion, we assume ti is fixed. If Pc{xq, Zi) is nonzero, then the 
state cannot reach Xi from xq by using admissible controls. The maximum 
value of Pc{xo, zi) over Dq x Di represents the worst scenario of reachability. 
In some applications, the relative ambiguity 

Pcixo,Xi) 

IkilU 

is used to measure the reachability. 

The definition of pc is consistent with the classic definition of controlla- 
bility for linear time-invariant systems. To be more specific, consider a linear 
system 

x^Ax + Bu (21) 

let u{-) E C be the space of continuous functions from [to,ti] to 3?"*; and let 
Dq^ Di^ 3?". If {A, B) is controllable, i.e. 



rank{ 



B AB A^B ■■■ A'^-^B 



n 



then for any xq and xi, there always exists a control input so that x{t) 
with x{0) — Xq reaches xi at t = ti. Therefore, pc{xo,xi) is always zero 
for arbitrary (xo,Xi); and pc = 0. On the other hand, if {A,B) is uncon- 
trollable, then under a change of coordinates an uncontrollable subsystem 
can be decoupled from the controllable part of the system. In the uncon- 
trollable subsystem, the states cannot be driven to close to each other by 
control inputs. Therefore, Pc{xq,xi) is unbounded for arbitrary states in the 
uncontrollable subspace. This implies pc — oo. To summarize, 

_ f if (A, B) is controllable 
^'^ 1 oo if (A, B) is uncontrollable 

Control has a cost. For weakly reachable systems, it takes relatively large 
control energy to reach a terminal state. The cost in reachability can be 
measured by the following quantity. Denote and \ \x\\x the metrics 

of the control input and the state, respective. 
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Definition 3 Given initial and final states, Xq and Xi, define 



W{xo,Xi) =min lim {\\u{t) - u*{t)\\u + '4'\\x{ti) -a;i||x) 



subject to 



X = f{x,u) 
x{to) = Xo, 
{x{-),u{-))eC 



(22) 



This definition has the following property. If xi can be reached from xq, 
then W{xo,xi) equals the minimum control 



If x{t) cannot reach Xi using admissible control, then W{xq,Xi) = oo. A 
large value of W{xo,Xi) implies higher control cost, thus weak reachability. 

Rcviark 7. Suppose the system is linear. Suppose Dq — {0} and Di — 
{x\ ||a;|| < e} for some small e > 0. Let W be the maximum cost in 
reachability under L^-norm, 



Then {W/eY equals the reciprocal of the smallest eigenvalue of the control- 
lability gramian. 

To justify Remark 7, consider 



If (A, B) is uncontrollable, we know W = oo. We also know that the smallest 
eigenvalue of the controllability gramian is zero. Therefore the claim holds 
true. Now suppose {A, B) is controllable. Define 



W{xq, Xi) = min | \u{t) - u*{t) \ \u 



{x,u) 



subject to 

X = f{x,u) 

x{to) = Xo, x{ti) = Xi 

{x{-),u{-))eC 



W = max W{0,xi) 



X — Ax + Bu 




to 



24 



The control cost to reach xi is defined by 

W{0,xi) = mm\\u{L)\\u 

subject to 

X = Ax + Bu 
x(to) = 0, x(ti) = xi 

Let P be the controllability gramian, then it is known [22] that the optimal 
cost satisfies 

ti 

\HtWdt 

to 

= xJP^^Xi 

If c"mm is the smallest eigenvalue of P, then 

1 



(W/ey = — max xTP xi 

e ||a;i||<<: 

1 



4.1 Example 

In the following, we compute the ambiguity in the control of a heat equation 
with Neumann boundary control. 

dw{r,t) d'^w{r,t) ^ 

dt dr'^ 
w{r, 0) = 0, < r < 27r (23) 

w(0,t) = 0, 0<t<tf 

Wrt{r,t)\r=27T = U{t) 

where w{r, t) G 3fJ is the state of the system, r G 3? is the space variable, and 
t is time. The control input is u. For a thermal problem, u represents the 
rate of the heat flux Wr- The initial state is assumed to be zero. 

For the purpose of computation, we discretize the problem at equally 
spaced nodes, 

= ro < ri < r2 < ■ ■ ■ < r^r = 27r 

Deflne 

xi{t) = w{ri,t), X2it) = w{r2,t), ■ ■ ■ ,a;7v(t) = w{rN,t) 
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Using central difference in space, fl23l) is approximated by the following con- 
trol system defined by ODEs 

Xi 
X2 

Xi 
Xn-1 

Xn 

V 

We understand that more sophisticated algorithms of solving the heat 
equation exit. We adopt this central difference method for simplicity in the 
illustration of control ambiguity. System (|2^ is linear and controllable. So, 
it is theoretically a reachable system. However, in reality the maximum 
temperature cannot exceed safety margin. Under such constraint, a control- 
lable linear system may not be reachable due to overshot. Let tf = 150 and 
K = 0.14. Suppose the target states are the following arches. A few of them 
is shown in Figure 14. 1[ 

Wf{r) = w{r,tf) = 74sin(r/2), < r < 277 

The goal is to compute the ambiguity in control from w{r, 0) = to Wf{r) 
with the magnitude of < A < 1.2 subject to the constraint 

w{r,t) < 2 

The norm in the finite dimensional state space is defined by 

N 

I |x| p = Ar ^ 

i=l 

which approximates the L^-norm in C[0,27r]. To compute the control ambi- 
guity in the range of < A < 1.2, we consider the following nodes in the 
magnitude of the target arch 

A = 0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.2 



X2 — 2xi 

Xi + X^- 2X2 



Xi—\ ~\~ Xi-\-\ ^Xi 



(24) 



xn-2 + a^AT — 2xn~ 
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Figure 8: The target state Wf{r) 



The corresponding dynamic optimization problem for the control ambi- 
guity is defined by 



TV 

2 



Pc(0,X/)2 = min Ar ^(xj - W/(rj)) 

subject to 

X = f{x,v) 

Xi{t)<2, i = l,2,---,iV 
x{0) = 

where f{x,v) is the dynamics defined in with the input v, and 



(25) 



Xf 



Wf{ri) Wf{r2) ■■■ Wf{rM) 



In the simulation, N is selected to be iV = 31. Problem (125|) is solved using 
Pseudospectral method at Legendre-Gauss-Lebato (LGL) nodes [T2 | [T^ [T5] . 
We use 15 LGL nodes in this example. Through computation, it is found 
that the system becomes increasingly unreachable due to the constraint when 
the value of A is bigger than 0.4. The relative ambiguity in control is shown 
in Figure 14.11 When the magnitude of the target state is 1.2, the relative 
ambiguity shows that the closest state that the system can reach has almost 
a 40% relative error. Therefore, the system is practically unreachable if the 
state is required to be bounded. 
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Figure 9: Magnitude of target state vs. relative control ambiguity 

5 Conclusion 

It is shown by numerous examples and definitions that computational dy- 
namic optimization is a promising tool of quantitatively analyzing control 
system properties. Using computational approaches, the concepts studied 
in this paper, including the ambiguity in estimation and control, input-to- 
output gain, and the cost in reachability, are applicable to a wide spectrum 
of apphcations. In addition, these concepts are defined and applied in a way 
so that one can take advantage of user knowledge or take into account sys- 
tem constraints. As a result, the properties of control systems are not only 
verified, but also measured quantitatively. While these concepts can be ap- 
plied to a wide spectrum of problems, some specific apphcations exemplified 
in this study include: strongly observable (detectable) or weakly observable 
(detectable) systems; improving observability by employing user knowledge; 
partial observability of networked complex systems; L^-gain of nonlinear con- 
trol systems; reachability in the presence of state constraints; and boundary 
control of partial differential equations. 

Similar to many nonlinear optimization problems, a main drawback of 
the approach is that a global optimization is, in general, not guaranteed for 
nonlinear systems. In addition, the problem of computational accuracy also 
poses many questions remain to be answered. 
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